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We investigate theoretically soliton excitations and dynamics of their formation in strongly cor- 
related systems of ultracold bosonic atoms in two and three dimensional optical lattices. We derive 
equations of nonlinear hydrodynamics in the regime of strong interactions and incommensurate fill- 
ings, when atoms can be treated as hard core bosons. When parameters change in one direction 
only we obtain Korteweg-de Vries type equation away from half-filling and modified KdV equation 
at half-filling. We apply this general analysis to a problem of the decay of the density step. We con- 
sider stability of one dimensional solutions to transverse fluctuations. Our results are also relevant 
for understanding nonequilibrium dynamics of lattice spin models. 

PACS numbers: 03. 75. Be, 32.80.Pj, 42.50.Vk 



INTRODUCTION 



Solitons are conspicuous manifestations of nonlinear interactions in a variety of physical systems (see e.g. [5ll l73l[). 
Originally introduced in hydrodynamics of classical fluids, they were later observed in a variety of other systems, 
including plasma physics, nonlinear optics, magnetism, dynamics of molecular systems. It is currently understood 
that formation of oscillatory zones and localized solitonic solutions is a common feature of many non-linear systems 
and does not depend on the exact integrability of the model. However the character of solitons is different for each 
system and understanding their properties remains a fundamental problem in physics and mathematics. 

In this paper we investigate theoretically the nature of solitons and dynamics of their formation in strongly correlated 
systems of ultracold bosonic atoms in optical lattices [H, 13, 42, Glj. Recently questions of far from equilibrium many- 



body dynamics took central stage in both theoretical and experimental study of ultracold atoms. What makes such 
systems particularly well suited for exploring quantum dynamics is their good isolation from the environment. Their 
characteristic energies and frequencies are of the order of kiloHertz, which is extremely convenient for experimental 
studies. It is also important that a wide array of experimental tools that allow to control system parameters in time 
and prepare far from equilibrium initial states have been developed. Recent experiments addressed such question as 
dynamics of fcrmions in optical lattices 0, [88| , observation of superexchange interactions using spin dynamics [oH , 



thermalization and relaxation in one-dimensional systems |34l . 1371 . 15011 . motion of impurity particles [531 , dynamics and 



adiabaticity in crossing classical and quantum phase transitions |78l . l8lj . Another important recent achievement is 



development of experimental tools for the in-situ imaging of individual atoms in optical lattices [ESZ! 45, 7^ 74. 87 1 



and low dimensional condensates 3^ 102 |. This technique allows unprecedented level of characterization of many- 



body states and should lead to deeper understanding of their out of equilibrium dynamics. One example is recent 
analysis of Bakr et al Q of the dynamics of defect creation in crossing from the SF to Mott state in two dimensional 
optical lattices. 

We start our analysis by deriving hydrodynamical approach to describe quantum dynamics of the lattice bosons. 
Hydrodynamical description has been applied to quantum many-body systems previously, including superfluids [48| . 
superconductors [3l)|, quantum Hall systems [93], and magnets [ssl]. The focus of most earlier analysis was on un- 
derstanding collective modes and universal features of linear response functions, which only required understanding 
linear hydrodynamics. When non-linear effects have been discussed for superfluid systems, it was primarily done for 
systems in the continuum with the the full Galilean symmetry. Our goal will be to include both nonlinearities and 
dispersion, since the competition of the two determines the nature and dynamics of solitons. 

Solitons in systems of ultracold atoms have been discussed previously in the regimes where semiclassical Gross- 
Pitaevskii equation can be applied either in uniform systems [3, [3 EM S, [z^ or systems with optical lattices 
0, S [H, 43, 47, 0, ll^. In this paper we will be interested in the regime of very strong interactions between 



atoms, the so-called hard core bosons regime [82|, |85|. In this case dynamics of atoms in a lattice can be described 
using anisotropic Heisenberg model [82l|. We demonstrate that in this regime the character of soliton excitations is 
very different and depends on both the filling factor and parameters of the Heisenberg model. Numerical analysis 
of solitons in Bose systems in optical lattices in the vicinity of the SF/Mott transition has been done recently by 
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Krutitsky et al. 55| . Our results can also be applied to study nonequilibrium spin dynamics of two component Bose 



mixtures in the Mott insulating regime 2l|, [52| and lattice spin systems in solid state physics. 



MODEL 

From lattice bosons to spin Hamiltonian 



Microscopic model describing ultracold bosonic atoms in an optical lattice is given by the Bose-Hubbard model 13, 

Si 

Hbh = -tb E ^I^J- + I E - 1) (1) 

(u) « 

Here b\ is a creation operator for bosons on site i, n.i = b\bi is the number of atoms on site i. We do not include the 
chemical potential term because in this paper we study dynamics and the operator of the total number of particles 
commutes with the Hamiltonian. It is sufficient to impose a certain number of particles at the initial time and then 
the total number of particles should not change during evolution. When there is inhomogeneous external potential 
we also need to add 



To keep the model more general we include nearest neighbor interactions 

■^Ext BH = 'HHub + ^ X] 

Such non-local interaction are relevant for atoms in higher Bloch bands [s^] and polar molecules in optical lattices [sg}. 
We consider a regime when the local repulsion U is large and the density of particles is incommensurate with the 
lattice. In this case strong number fluctuations are suppressed even in the superfluid state and we can limit the Hilbert 
space to only two possible occupation numbers |rto ^ 1) ^^nd |rto) per site. It is convenient to represent these states as 
spin states. State jno — 1)^ corresponds to | J,)^, and state |no)i corresponds to | In this limit Hamiltonian ([3]) is 
equivalent to the anisotropic Heisenberg model 

Hah = - J± E + -J^H ^'^l (4) 

Here tr" are Pauli matrices, 2J± = t^no, and Jz — —V- 

Hamiltonian ^ also appears as an effective description of spin dynamics in the Mott state of two component Bose 
mixtures at filling factor n = 1 [2l|, . 

Semiclassical equations of motion for lattice bosons 

In this section we discuss how one can obtain semiclassical description of dynamics of ^ using either variational 
Gutzwiller wavefunctions or linearized equations of motion, which in this case are equivalent to lattice Landau-Lifshitz 
equations. To simplify the derivation we assume that parameters of the system change in one direction only. We 
emphasize that our focus is on two and three dimensional systems. Restriction to having variations of parameters in 
only one direction is, firstly, for notational simplicity (extension to higher dimensions is straightforward) and, secondly, 
because we will be concerned with problems where initial state has been prepared to have parameters changing along 
one of the coordinates. We discuss effects of fluctuations in transverse directions in subsequent sections. 

Strictly one dimensional systems are special and mean-field approaches do not apply to them even in equilibrium. 
However s pec ial analytical approaches are available for one dimensional systems, including fermionization and Bcthe 



ansatz[25|, |29|, |89|. Also powerful numerical methods based on DMRG 86] and Matrix Product States[l7| allow to 
study dynamics of one dimensional systems in great details. On the other hand, nonequilibrium dynamics of higher 
dimensional systems remains largely unexplored. This is the main motivation for the current paper. Interestingly, 
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recent work by Lancaster and Mitra 57| showed that semiclassical analysis of Landau-Lifshitz equations for one 
dimensional spin chains give results consistent with exact calculations. Hence our results may also be relevant for one 
dimensional spin chains. 

To obtain semiclassical dynamics we consider time-dependent variational wavefunctions 



= n 



. 0^{t) 



COS ■ 



0^{t) 



(5) 



Expectation values of the original boson operators are 



im) = no + - {cosOi - 1) 



To project Schrodinger equation into wavefunction 1^ we define the Lagrangian 4^, 41 1 



(6) 



j:1 



— J± sin^i sin6'j cos{(pi — ^Pj) — Jz cos6'i cos( 



and write equations of motion 



d 5L _^ 
dt 6qi 6qi 



Here qi corresponds to both (pi and 9i. We find 

ifi sin 6*,, = 



= -4 Jj_ cos6lj (sinfi^i+i cos(i^i - ipi+i) + sin&i^i cos((/3j - (fi-i)) + 
+ AJz sin 9i (cos(?i+i + cos0i_i) 

= -4 Jj_ (sin6'i+i sin((pi - 'Pi+i) + sin6'i_i sin((pi - 'Pi-i)) 



(7) 



(8) 



The first equation is effectively the Josephson relation: time derivative of the phase is equal to the chemical potential 
which depends on the values of 9 and Lp. The second equation is charge conservation. 

One can give an alternative physical interpretation to equations (|S]) . We write equations of motion for spin operators 



(9) 



And we have analogous equations for a, 
tation values 



To obtain semiclassical dynamics we replace operators by their expec- 



dt 



(10) 



These are familiar Landau-Lifshitz equations. If we use wavefunction to calculate {a^) — cosOi and {erf) 



: sm f,e 



we recognize that Landau-Lifshitz equations are equivalent to (|H]). 



Dyna mics of the Bose-Hubbard model has been studied using Gutzwiller variational wavefunctions in [18|, |39| . 
H, [lOll. In 0, [ttI this approach was used to describe current decay in the strongly interacting regime of bosons. 
Theoretical predictions were in quantitative agreement with subsequent experimental results by Mun et al[67[. 



SEMICLASSICAL DYNAMICS IN THE CONTINUUM LIMIT 



Long wavelength expansion 



It is convenient to introduce slow variables in space, X = hx, and time, T = ht, where h is the lattice constant. We 
are looking at dynamics of fluctuations that are slow on the scale of the lattice constant. So /i is a small parameter 
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in which we will expand. We introduce 



and obtain 



fjL = cos 

a(X,T) = hip{X,t). (11) 



1 U^U^ 

£ = -(TT - 2 J± (1 - fi'^) COS ax - 2 J^fi'^ + h'^ J± I-Jz2L cos ax + 

2 1 — /i^ 

h^J^{l-y?) \ \axxx sinax + 7 coscrx ) - Jz^J-^J-xx + (12) 



Hydrodynamics 

If we keep only the lowest order terms in h in (jl2p . we obtain the hydrodynamic part of the lagrangian 

C-Hydr = -CTTM - 2 J_l (1 - COSCTX - 2 Jz^^ (13) 

It is convenient to define 

k{X,T) = ax{X,T) 



(14) 



The new variable is proportional to the phase gradient, k ~ V(^. 

Equations of motion obtained from the lagrangian (|13p have a standard hydrodynamic form 



kx = 8 J± ^J, sinfc kx — (8 Jl cosfc — 8 J^) /ix 
Ht — — 4 J_L (1 — /i^) cos fc fcx + 8 J± n sin fc /xx 



(15) 



Linearized equations of motion. Stable and unstable regimes. 

Let us consider a superfluid state with a uniform density and, possibly, finite phase winding. When fc 7^ 0, this is 
current carrying state with / = 4Jl(1 — ^q) sinfc. 

Frequencies of linearized excitations are given by the eigenvalues of the matrix 



. _ , 8 Ji/i sinfc — 8 Jj^ cos fc + 8 . 

~ I -4:Jj_{l- ^i^) coak 8Jj_^iamk ' ^ ' 



We have for the eigenvalues of A 



Ai,2 = 8Jj_fi sin k ± \/4 J_L (1 - /i2) cos fc (8 J± cos fc - 8 J^) (17) 

When Jj^ > Jz and fc is small, both eigenvalues of ([T6l) are real (when Jz — this is true for all fc ). This is the 
hyperbolic regime, which will be the main focus of our paper. 

When < cosfc < Jz/J±, eigenvalues of appear as a complex conjugate pair. This is the elliptic regime, 
which corresponds to the unstable state of the system. In this regime small fluctuations of the plane wave type 

fc(X,T) = fco + SkiX,T) , ii{X,T) ^ /io + SiiiX,T) 
5k{X,T) - 5ke^i^+"'^i^'^ , 5ii{X,T) - s^^(,^qX+iu(q)T 

grow exponentially in time. Existence of this unstable regime is known as the dynamical instability [i, 94|. It was 
observed experimentally for atoms in optical lattices 23|, |67| . Exponential growth of small modulations predicted by 
equations ([T5|) is only valid for short times. Dynamics of the unstable regime beyond the short time limit can be 
analyzed using mathematical methods from the theory of elliptic equations. In this paper we only address the stable 
hyperbolic regime. 
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When the initial state does not carry a current, i.e. there is no phase winding, 

ipiX,0) = (18) 
To obtain further insight into the hnearized system we set 

^,iX,T) = /io+/7(X,T) 

We can now rewrite equation ([T5|) in terms of variables p{X, T) and k{X, T), which describe small deviations from 
the equilibrium state 

kT = -8{J±-J,)px , PT = -^J±{l-Po)kx (19) 
System gives the following equation for p{X, T) 

PTT = ■i2J^{l^ pI){J^- J,)pxx (20) 

We find the familiar wave equation, which describes propagation of the initial perturbation p{X, 0) with a small 
amplitude. Equations (jl9p and (|20p show that during the dynamical evolution of the perturbation, the superfluid 
velocity k{X,T) is of the order of p{X,T), provided that this is true in the initial state. This is the regime that will 
be the focus of our paper. 

Nonlinearities and appearance of singularities 

We now include nonlinear terms in the analysis of equations of motion. In the simplest case = we can define 

— \pi arcsin/i — fc = — \f2Q — k 

= V2 arcsin^ + k = n/Vo, - V2e + k (21) 
And from the Lagrangian (|13p we obtain equations of motion 

= (S Jj_ sin^^ sin^^ + 4^2 Jj_ cos cos 
4 = {sJi_ sin 4^ sin - 4V2Ji_ cos cos 



r.2 ^1 



2 I 'X 

(22) 



r.2 ^1 



—J 'X 

Equations (j22p are written in terms of the Riemann invariants, which separate the system (jisp into the left- and 
right-moving parts. This representation is most convenient in the analysis of Hydrodynamic Type systems. System 
of equations (1221) admits two natural reductions = const or ~ const, which describe separate propagation of the 
left- and right-moving excitations. 

When Jz ^ 0, expressions for Riemann invariants are more cumbersome 



r 



/- . fK ■ Vcosfc 

= V 2 arcsm /i — v 2 arcsm /ig T / — 1^^^=^^^= d-k 

Jo \/ cos fc — Jz/ J\ 



-y/cosfc — Jz/ J± 



The corresponding diagonal system of the equations of motion has a character close to (|22|) in the hyperbolic regime. 
Taking in the account that the functions p{X, T) and k{X, T) have the same order in our approach we can write 

A1.2 = ± x/32 (1 - ^g) (Ji - Jz) + SJ^fiok T -1^ v/32 J±(J±- J.) p + 

V l^'^o 



(23) 

In the same way 



+ 8J^pfc T V2(l-A*^)^± T ^8J^{J^-Jz) ,' p^ 
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T=To 


< T < T(, ^~~X 


T = 0/^\ 













FIG. 1: Formation of the breaking point for equations (|26p . H27|l in the case /xq > 0. and describe left and right movers 
respectively. Dynamics is shown in moving frames of references. In both cases the rear edge of the wave steepens and develops 
a singularity. 



„i^2 _ V2 , V2^io p2 ^2(1 + 2^ll) p3 



To understand the role of non-linear effects we expand the corresponding equations of motion up to second order 
terms in deviations from the uniform state. We obtain the following general form of the equations of motion 

A A (25) 



X 



Equations \2b\ describe coupled evolution of the right and left moving parts. To get further insight into dynamics we 
make another simplification. In the problems that we consider the left and right moving parts overlap at short times, 
but separate after a finite time. The main effects of non-linearities appear at long times. Thus when discussing effects 
of non-linearities it is sufficient to consider separately the left- and right-moving parts of the solution. So when we 
discuss the dynamics of we can set = const and vice versa. 

After we make the Galilean transformation for the left and right propagating parts we obtain 



^ - ^^^oVJZ{JZ^J^)r^r], (26) 



4 = &tJ^o^J^{J^-Jz) r'^rl (27) 

Equations (|26|) and ([27]) are known as the Hopf equations describing " simple waves" . Their solutions are given by 
the implicit formula 



The most important feature of these solutions is that they exist only up to a finite time Tq, which depends on the 
initial conditions. All nontrivial solutions become singular after some finite time. Physically this corresponds to 
formation of the breaking point, which we show in Fig. [T] This can be understood as a result of regions of different 
densities moving with different velocities. 

Formation of the singularity is not restricted to the truncated equations of motion ([25]). This is a feature of the 
general non-linear dynamics of the equations of motion (j22p . Generally system of equations (jl5l) can be reduced to 
a linear problem using the so-called hodograph transformation. Then solutions of (|15p can be described in terms of 
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perturbations moving along the characteristic hnes dX/dT = Ai^2(?'^, J'^)- Characteristics of the nonhnear system 
depend on the variables (r^,r^) and unique solutions of (jlSp exist only up to a finite time Tq. At later times solution 
becomes multi- valued. Special solutions of (|15p given by relations r^(X, T) — const or r^{X,T) = const describe 
perturbations moving along one of the characteristic lines. In this case the second variable (r^ or r^) satisfies a 
nonlinear first order equation, which is (locally) equivalent to the nonlinear Hopf equation. 

For times approaching Tq solutions of (^5)) . ([27]) are close to developing a breaking point and have high gradients. 
In this regime neglecting higher order gradients in the Lagrangian (|12[) is no longer justified. In the next section 
we will see that taking dispersion into account suppresses singularities in the solutions and gives rise to short-period 
oscillations. 

General analysis of how dispersion leads to the formation of oscillatory zones in our system is rather complicated. 
In the most generic case, one can not use expansion (1121) to describe the oscillatory zone formation. The period of 
oscillations arising for T > Tq is of the order of h, so all higher dispersive corrections are of the same order. Accurate 
description of the transition from the " slowly- modulated" to the rapidly modulated regimes can only be done with 
the use of the original lattice system (|8]). However, there are certain special cases, in which the use of the continuum 
model (fT2| is justified. Fortunately these cases are interesting from the experimental point of view. They will be the 
subject of our discussion. 

NONLINEAR WAVES IN GENERIC CASE 
Connection to Korteweg-de Vries equation 

When discussing dispersive terms for p and k in the equations of motion, it is sufficient to keep them only in the 
linear order in deviations from the uniform state. Dispersive terms come with additional factors of h and are already 
small. Hence in the Lagrangian (jl2p dispersive terms need to be considered only up to quadratic terms in p or ax- 
Modulo total derivatives with respect to X we can write 

r J.2 T /^O 2 , 1,2 T 2 

i-Disp — n J± ^ px + ri Jz Px 

^ Mo 

The resulting equations of motion are 

kr ^ 8J±p sin kkx - {8 J± cos fc - 8 J^) px + 4:h'^ J± pxxx + 'ih'^ Jz Pxxx 

(28) 

Pt — — 4 Jj^ (l — ^^) cos k kx + 8 J_L /i sin/c px — \ h"^ J_l{1 — Pq) kxxx 
Using variables r^{k,p), r^{k,p) we obtain 

4 = ^/jlUr^Jl) (x/32(l-/z^) - 6por^ + 2por^) r], - 

-^h^ \/J±{J± - Jz) {rj,xx - rjcxx) " (29) 

-V2h^ [j±j^^+ Jz) v/W^y^5: Hxx + r%xx) 

4 = VJ±iJ±-Jz) (- v/32(l - Mo) - 2por^ + Qpor^) 4" 

-J^h^ VMJi- - 'h) VI -Mo irlxx - r]cxx) + (30) 

As in our earlier discussion we consider separately the left and right moving parts, i.e. we take either = const or 
= const. After we included effects of dispersion such reductions are no longer exact. However, in cases of interest. 



- — /i^ J_L (1 - pI) a\x 
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interaction between and gives rise only to small rapid oscillations. Such oscillations are expected to be much 



smaller than the structures that we discuss (see e.g. 5j|) and we will neglect them in this paper. We also perform 



Galilean transformations for the two parts and obtain in the moving coordinate systems 



4 = 6^.0 ^JAJ^-J.) r\ - V2h' y'l - ^4 (j± - - \ J.) r\^x (32) 

Note that equations (PT|) . (|32p transform into each other if we change X ^ —X. Equivalence of the two equations 
for fixed values of Jl, and /io represents an evident corollary of the symmetry X — —X of the original system. 
It is not difficult to see that equations (l3Tt . (|32|) represent the KdV-equation provided that 



A^o ^ , J± ( - - ) - - ^ 



.6 6 

Depending on the values of parameters Jl, J^, and /io, equations pip and p2p are equivalent to one of the following 
two equations 



{/t + 6C/t/x - C/xxx = (33) 

C/t + 6 [/ J7x + Uxxx = (34) 

after an appropriate rescaling of coordinates (X, T) and functions . These two equations are equivalent to each 
other if we admit the inversion — >■ — r% T — T. Howev er, t his transformation leads to very different physical 
interpretation of solutions for a fixed /ig, as we discuss below. [lOSj 

KdV type equations and allow solitonic solutions, which are long lived nonlinear excitations in the system. 
The velocity of a soliton is proportional to its amplitude, so larger solitons move faster than the smaller ones. The 
asymptotic form of an iV-soliton solution for T — >■ oo for equations ([55)1 and ([M)) can be represented as shown at Fig. 
[2j In the Appendix □ we briefly review how one can verify the existence of solitonic excitations in the KdV equation 
using connection to the linear Schroedinger equation. We also point out that in general, solutions of ([33|) and (p4|) 
include not only the soliton part but also " wave trains" . The soliton part and the " wave train" parts separate from 
each other at long times (Fig. [3]). The soliton part of the solution remains unchanged for all T > while the wave - 
train part "dissolves" as T — >■ cxo ([s^). From our point of view, solitons of ([55]) and (|34l) represent the most interesting 
part of the solution and we focus on them in this paper. 



Discussion of solitonic excitations 



Solitons in KdV equations have been studied in detail during the last few decades. In this paper we take previously 
known mathematical results and discuss their physical implications for our specific system. While we provide a brief 
summary of the mathematical methods used in analyzing soliton excitations in the Appendix, we refer readers to the 



books PJ, [TlJ, [T^l for a more detailed discussion of general mathematical aspects of the KdV equation. 

The character of solitonic solutions of KdV type equations and depends on parameters. In particular 
depending on the ratio of Jz/J_l and the density, isolated solitons can appear either as particle-like or hole-like 
excitations. In this subsection we only provide a summary of the results. More details can be found in the Appendix. 

In the discussion below we only consider the case /io > 0, which corresponds to the density above half-filling 
(n) > 1/2. Equations (|3ip and (|32p have a symmetry /zo — > — /xo, Ti — ^ — r^. This symmetry originates from the 
particle-hole symmetry of the initial system, which relates states below and above half-filling ^ : B ^ -n — B and 
<y9 — > — In our discussion this symmetry allows to relate solitonic excitations below and above half-filling. For 
/io < solitons are "mirror images" of the /^o > case. For example, if we find particle-like solitons above half- filling, 
we should have hole-like solitons below half-filling (/^o — ^ — /^o) for the same values of J. Let us represent here also 
the form of the "hole- like" and the "particle- like" solitons in the original variables (k, p) (see Fig. [4]). 

We also remind the readers that we only need to consider states that are stable against dynamical modulations, i.e. 
J_L > Jz- 




FIG. 2: The asymptotic form (T — >■ oo) of the A^'-soliton solutions for equations H33[) and ((34} respectively (Vi > V2 > V3). 



Solitons for Jz > J±/7 and /lo > 0. 



In this case both equations ([5T|) . ([5^ reduce to equation (IMl) after rescahng the variables and, if necessary, per- 
forming the transformation X — > —X. There should be no solitons when U{X) < 0. We are guaranteed to find 
solitonic excitations vifhen 

/+00 
U{X) dX > (35) 
-co 

When U(X) > the soliton part represents the main part of the solution. So we find particle-like solitons in this 
situation. 



Solitons for < J±/7 and < fio < \J Z j 



) 

Now both equations (|3T|). ([32|) reduce to equation (|33|) after rescaling the variables and doing the transformation 
X —X in equation ([5^ . This equation does not have any solitons when U{X) > 0. It has guaranteed solitonic 
solutions when 

+ C30 

U{X) dX < (36) 

-00 

In the case with U{X) < the soliton part represents the main part of the solution. Hence in terms of the original 
density, wc find the hole-type solitons in this case. 



Solitons for Jz < Jx/7 and < fio < 1 



Both the equations pip . p2p reduce to equation ([M)) in this case. Thus we find particle-like solitons in terms of 
the original density. 





FIG. 4: The form of a soliton solution in the (k, p)-variables for the case of the "hole-Uke" soUton and the "particle- like" soliton 
respectively. 



We can now summarize results of this subsection. When > J±/7 we find that above half- filling there are only 
particle- like solitonic excitations. When < J±/7 and above half-filling we find that we have either hole-like (closer 
to half-filling) or particle-like solitons (closer to filling factor one). 



Self-consistency of the long wavelength expansion 



Before concluding this section we would like to verify that our solutions do not take us outside the region of 
applicability of Lagrangian (jl2[) . which was obtained using long wavelength expansion. When we consider dynamics 
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starting from a state with small smooth deviations from a uniform density, approximate Lagrangian (|T2)) can be 
certainly used at the initial stages of the evolution. However, at final (asymptotic) stages of the evolution, the 
solution may be sufficiently different from the initial state. Let us consider specifically the soliton part of asymptotic 
solutions. In soliton solutions both the nonlinear and dispersive parts are important and the interplay of the two gives 
rise to a stable soliton. One of the important properties of the KdV equation is that the amplitude of solitons is of 
the same order as initial deviations from the uniform density. Equations (l5Tt - ([5^ were obtained assuming small 
deviations of the initial density from the uniform value /^o • These small deviations set the scale for the amplitude of 
resulting solitons. In solitons there is a direct relation between the amplitude and the width (the width increases as 
the amplitude goes to zero). Hence in the limit that we discuss, the dispersion part of our soliton solutions should be 
small, and our approximation of neglecting higher dispersion corrections should be justified even for the final stages 
of the evolution. For example, when solution can be written as the " quasiclassical solution", in which the soliton part 
represents the main contribution to the solution, higher dispersive and nonlinear terms should have very weak effect 
on the soliton. 

Similar considerations are applicable for the "wave-train" part of the solutions. However, the "wave-train" part 
dissolves in the limit T — >■ oo and we expect that it will be more challenging to observe it in experiments. 



NONLINEAR WAVES IN SPECIAL CASES 



Half-filling. Solitons of the modified Korteweg-de Vries equation 



When the particle density is 1/2, the system of hard core bosons has a full particle-hole symmetry. Eigenvalues 
Ai,2 of the linearized system (fT9|) have the largest possible magnitude 



Ai,2 = ±4V2J±(J±- J.) 

which corresponds to the largest possible velocity of linear waves. In the case of dynamics starting from some initial 
state, this should provide fastest spatial separation of the left- and right-moving parts of the perturbation. In this case 
/iQ = 0, so corrections to Xi^2, which are linear in p and k, vanish and we need to use quadratic terms in the expansion 
([23|) . In our discussion below we keep linear terms, in order to accommodate small /io ^ 0. Using approximation ((24|) 
we can write 



Ai ~ a/ J_l(J_l — Jz 



4^2 - 6nor^ + 2nor^ 



2^/2J_ 



z 1 2 

r r 



Ai ~ ^ J±{J± — Jz 



-4\/2 - 2fior^ + 6nor^ 



J-L ~^ Jz / 1^2 , "^J-i- Jz ( 2\2 'J ^ '^z 

[r ) -i p [r j - 



2V2J± 



2V2J± 



2V2Jj_ 



From the last two equations we determine how propagation of the left- and right- moving parts, pi p - p2p . is modified 
by the higher order terms. Within the assumptions of spatial separation of the left- and right-moving parts, which 
we used in the earlier discussion, and using appropriate moving frames of reference we find 



4 = ^J^{J^-Jz) (e^or^ + -X ~ ^^(i^) - ^J^)-xxx (38) 

In writing the last equations we omitted higher order corrections in fiQ. When J± ^ 7Jz- equation (|37|) can be written 
in the canonical form 



Ut + {aU + Q U^) Ux ± Uxxx = 



(39) 
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after a scaling transformation. Parameter a that we introduced here is proportional to the deviation from half-filling, 
a /io, and we assume it to be small. 

Equation (|39p is called the modified Korteweg - de Vries (mKdV) equation and represents an integrable system as 
well as the KdV equation (see |93|). Let us note also that the mKdV equation is connected with the KdV equation 
by the Miura transformation ([66|) which was the first observation of the integrability properties of the KdV equation 



itself (see [71 



There is a wider variety of soliton excitations that one can construct in the mKdV problem. At a fixed value 
of the chemical potential one can find both particle-like and hole-like solitons moving in the same direction. This 
should be contrasted to the situation away from half-filling, which we discussed in the previous section, where at 
a given chemical potential and direction of propagation we had either particle or hole like solitons, but never both 
simultaneously. For the mKdV case we also find soliton excitations which look like particle on a pedestal (or hole 
on a pedestal). We provide a detailed discussion of solitons in the mKdV problem and their manifestations for our 
system in the Appendix. 

Close to integer filling. Nonlinear Schroedinger equation 

When the system is close to integer filling /io = ±1- In this case characteristic velocities of the linearized system 
(|19p coincide. Hence we can no longer assume separation of the left- and right-moving parts. Examining system 
((29)l - (|30p we find that dispersive corrections also have singularities in variables (p, fc) . 

To avoid these difficulties we return to variables (6*, cr), which we used before, and consider the Lagrangian density 

£ = ^ or cos6' — 2 Jl sin^ 6* cos ax — 2 cos^ 6* — h? J±^ sin 6* (sin6')j^j^ cosux — 

(40) 

— h'^J±_ sin^ (^(Txxx sinax + jcrjcx cosctx) — h'^ Jz cos9 {cos9)^-^ + ©(/i'') 

in the limit — > 0, vr. If we keep only quadratic terms in {sm6,a) in the dispersive part of the Lagrangian, we can 
write the corresponding equations of motion as (in the limit fia = ±1) 

ar sin 9 + 8 J± sin 6 cos 9 cos ax — 8 Jz cos 9 sin 9 + 4h^ J± fiQ (sin 6*)^^ = 

(41) 

9t ^ 8 J± (sin 9)j^ sin ax — 4 Jl sin 9 (sin ax)x — 



If we keep only the lowest order cubic terms in the nonlinear part of (1411) , we can rewrite this equation for small 
{9, a) as 

ar sin ^ -I- 8 /zq ( J_l — Jz) sin 9 — A fio ( Jj_ — Jz) sin^ 9 — 4 fia Jj_ sin 9 aj^ +4 h"^ fio J± (sin 9)xx — 

9t — 8J± (sin 9)j^ ax ~ 4 Jj^ sin 9 axx = 

(Mo = ±1). 

It is not difficult to verify that the system above can be written in the form of the defocusing nonlinear Shrodinger 
equation 



ihiPr = 8 (Ji - J,) V - 4 (Ji - Jz) IVIV + 4h^J^ iPxx (42) 



for the function 



= sin 61 e'^"'"/'' = sin 9 e*'""'^ 

Equation p2)) describes an integrable system ([99]), which was solved by V.E. Zakharov and A.B. Shabat by the 
inverse scattering method. System (|42|) admits an exact description of the evolution starting from any initial state. 
However, nonlinear Shrodinger equation does not have soliton solutions in the defocusing case. Defocusing nature of 
equation (j42p demonstrates stable behavior of the system with respect to initial perturbations. In this case asymptotic 
behavior of solutions of (|42l) should only include wave-trains which "dissolve" for T — > oo (js^). |l04j 
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FIG. 5: The step-like initial conditions for the function p[X, T) with the assumption tp(X, 0) = 0. 



Special filling factor 

We now comment on the special point of our system at 

Mo = ^Ji_ - IJz I V7{J± - Jz) (43) 
for the case < < J_l/7. To get equations ([55)1 - from (PT|) - ([5^ we need to make scaling transformation 

This transformation is singular at the special point (j43|) . As a corollary, the width of solitons (and the period of 
oscillations in the "wave-train" part) become small in X-space w.r.t. another parameter 

Mo - VJ± - / V7iJ± - Jz) 

Higher dispersive terms become important in this limit and Lagrangian density (|12p can no longer be used. As we 
discussed earlier, dynamics is more complicated near this special point, and one should use original lattice system (|5]) 
to discuss dynamics. In general we expect here oscillation zones with rather short period of oscillations. 

DECAY OF THE DENSITY STEP 

We now apply our general arguments to understand dynamical evolution starting from a specific initial state. We 
assume that at T = we have a smooth step-like change in the density without any initial current. Experimentally 
such initial configuration can be created using a smooth step in the external potential that is suddenly removed. This 
initial state is of the form given by equation (|18l) . It is shown schematically in Fig. [SJ In this section we only consider 
the situation when the system is not close to any special points. Density step decay for systems close to half-filling is 
discussed in the Appendix. 

Since we rely on the long wavelength expansion, we assume that function p {X, T = 0) is a slow function of the 
spatial coordinate. 

The main terms in the long wavelength expansion of dynamics are given by the wave equation (jl9p . The wave 
equation predicts that after a short time the step-like initial state should turn into a two-step solution, with two 
steps propagating in the opposite directions (see Fig. [6]). When the two steps separate from each other, they can be 
analyzed independently. The left- and right-moving edges of the solution correspond to ri and r2 . Proceeding to the 
next order in h, we find that they are described by equations ([3T1) and ([32|) respectively. After rescaling of coordinates 
{X,T) and functions and themselves, this dynamics is given either by equation (1551) or where the choice 
depends on the values of {J±, Jz, and mo)- 
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FIG. 6: The evolution of initial distribution with k(X, 0) = and the step-like p(X, 0) in the approximation of system (|19p . 




FIG. 7: The increasing of the steepness of solution r^(X) on the left-moving edge and the decreasing of the steepness of solution 
r'^(X) on the right-moving edge in the hydrodynamic approximation. The data are sketched in the original coordinate system. 

First of all, we need to understand whether hydrodynamic solutions for r\ and r2 break down and develop a 
singularity. For /io > and the initial density profile shown in Fig. [5] the steepness of function r^iX) should increase 
with time while the steepness of function r^iX) should decrease with time (see Fig. [7]). This follows from simple 
hydrodynamic analysis following equations (|26l) and (|27)) . This means that the steepness of solutions p(^), k{X) 
will increase on the left-moving edge of Fig. [6] and decrease on the right-moving edge. (The situation changes to the 
opposite for the inverse step initial state.) So in this case no dispersive corrections are needed for r^(X). Function 
r^(X) should remain smooth for all T > in the hydrodynamic approximation. On the other hand, function r^(X) 
develops a breaking point in the hydrodynamic approximation. Thus we need to consider equation (j3ip taking into 
account dispersive corrections. As we discussed before dispersive corrections should give rise the oscillation zone, 
which we expect to grow linearly with time. The form of oscillations should be different for equations (|33p and ([34]) 
due to different signs of dispersion in these systems (see Fig. |8]-[9|). 

To understand the oscillation zone that arises following breaking of the hydrodynamic solution we need to analyze 
dynamics of the KdV equation with step like initial conditions. This problem was addressed by A.V. Gurevich and 
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FIG. 8: The development of oscillation zone from the step- like initial data for the case of equation p3[) . 




FIG. 9: The development of oscillation zone from the step- like initial data for the case of equation (|34|) . 

L.P. Pitaevskii (fsil, 'ssj) using the Whitham theory of slow modulations. We will now summarize their key results 
pointing out their implications for our system. 

Gurevich and Pitaevskii considered a general problem of slowly modulated one-phase solution of the KdV equation. 
One-phase solution is a periodic running wave solution that provides a generalization of the one-soliton solutions of 
the KdV equation 

U{X,T) = <i>{KX +ujT + eo, K,A,n) 

One-phase solution depends on three parameters (K,A,n). Functions ^{0, K,A,n) should be 27r-periodic in 9, so 
parameter k plays the role of the wave number for nonlinear running waves. Parameter A plays the role of the 



<i>(0,K,A,n) 
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FIG. 10: The general form of the function $(0, k, A, n) representing the one-phase solution of the KdV - equation. 

amplitude of the periodic solution, while parameter n =< <i> > is the value of $ averaged over one period (see Fig. 
One-phase solutions of KdV can be written in the form 



^(kX + ojT,K,A,n) = dn 



where s is the modulus of the Jacobi elliptic function dn(M, s), < s < 1. The values {K,uj,n) can be expressed in 
terms of the parameters {A, s, 7) in the following way 

f A 47r / A / A 



A 



Us 



1/2 

2, {X-VT),s 



+ 7 



AE(s) 
s^K(s) 

where K{s) and E{s) are the elliptic integrals of the first and the second kind respectively. 
We can also write 



HO,A,s,j) = ^dn^(^^0,.) +7 

as normalization of function $(0, k, A, n). 

The one-soliton solutions of KdV can be considered as the limiting case of the one-phase solutions in the large-period 
limit K — 0. Traditionally the asymptotes $(0) — 0, 6* — ±cx) is assumed for the soliton solutions of KdV, so the 
amplitude parameter A remains the one parameter of a one-soliton solution. 

In Whitham's approach parameters (k. A, n) become slow functions of x and t 



K = k{X,T) , A = A{X,T) , n = n{X,T) 

so that functions k{X,T), A(X,T), n{X,T) satisfy a nontrivial system of quasilinear equations in partial derivatives 
(the so-called Whitham's system). Whitham's system describes evolution of initial parameters k{X,0), A{X,0), 
n{X, 0) of oscillating solutions, such that development of oscillations can be calculated in this case. 
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Gurevich and Pitaevskii showed that in the KdV equation with a step hke initial conditions, the small oscillation 
zone, that arises near the breaking point of the hydrodynamic solution, can be described by the self-similar solutions 
characterized by only one variable, I = X/T. 

In more details, the asymptotic (T — >■ oo) form of oscillations can be described by the modulated one-phase solutions 
of KdV with parameters k{X, T), A{X, T), n{X, T) of the form 

k(X,T) = k{X/T) , A{X,T) = A{X/T) , n{X,T) = n{X/T) 
The oscillation zone is located in the interval 



< X/T < 1+ 

in this asymptotic regime. 



According to [3^ - [33' the amplitude of oscillations A{X,T) becomes zero at the "trailing edge" of the oscillation 
zone (the right edge in Fig. |8]and the left edge in Fig. |9]). The wave number of nonlinear oscillations k{X, T) becomes 
zero at the "leading edge" of the oscillation zone (the left edge in Fig. [5] and the right edge in Fig. (5]). 

We can see that the "trailing edge" of oscillation zone can be considered as a source of oscillations with small 
amplitude, which develop into solitons in the limit T — )■ oo. The "leading edge" of the oscillation zone can be 
considered as a source of free solitons since we have k — >■ on this edge and the distance between solitons tends to 
infinity for T — > oo. 

We point out that it is also possible to analyze the problem above in terms of the "pure" soliton picture ( [sol. [60| ) . 
Approach used in ([1^, 11^) is also a classical part of the soliton theory. 

General problem of the decay of different initial configurations in the theory of small-dispersion KdV-equation 
represents a big branch of the soliton theory. While we do not discuss other problems here, we expect that many of 
the known mathematical results will be relevant for different experiments with ultracold atoms. We also point out 
that our methodology for identifying the character of solitons (particle- or hole-like) was based on considering the 
function U{X,T), which describes Riemann invariants r^{X,T) or r^{X,T). It is more natural to classify solitons 
based on the density. Relations between A^'^^ {X,T) and the more physical variables of the density, p{X,T), and 
the phase gradient, k{X,T), are given in equation (1^^ . We find that the density always follows the behavior of 
A^'^^ {X,T). Hence our classification of the hole-type and the particle-type solitons in terms of the density coincides 
with that given in terms of the function U{X,T). 

Before concluding this section we would like to point out that whether step- like conditions shown in Fig. [5] should 
be considered as a source of hole-like or particle-like solitons in the solutions p{X, T), k{X, T) depends on the relation 
between parameters { J±, Jz, fJ-o)- In general, we expect that larger values of Jz and fiQ suppress the appearance of 
hole-type solitons and favor solitons of the particle type. On the opposite side, smaller values of Jz and /io allow 
solitons of the hole type and suppress solitons of particle type. 



TWO-DIMENSIONAL EFFECTS. 



In this section we discuss the role of transverse directions. We consider a question of whether one dimensional 
profiles, that we discussed so far, are stable against "weak" modulation in the transverse direction. 

For a Z)-dimensional lattice we need to change the long wavelength Lagrangian density (|40|) to a more general 
expression 



1 ^ 

£ = - (Tt cos — 2 J± sin^ 9 cos cr^i — 2 Jz D cos^ 



i=l 



D 



J± sm0 ^ {sinff) ^i-^i cosaxi — h? Jz cos0 ^ (cos( 



J± sin^ 



D 

E 

1=1 



-(Tx-xixi sincTjfi 



ia^,^. cos ax.) + 0{h^) 
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FIG. 11: Schematic sketch of a one-dimensional sohton string modulated in the Y-direction. 



or 



C 



2 J_L (1 - fi'^) ^ cosaxi - 2J,_Dfi'^ 



1 - /i2 ^ 



D 

^ fix- COSCTxi 



D 

1=1 



D 



4=1 



J± (1 ^ A*^) X! ( q ^xixix- smax- + - o\ix- cos ax- 



in the coordinates (/x,cr). 

We separate the hydrodynamic and dispersive parts of the Lagrangian and repeat considerations used in the previous 
sections. Analysis of the dynamical system is more complicated for D > 1 and we will not explore all of its richness. 
We only address a question whether one dimensional solitons that we discussed so far are stable with respect to 
formation of a two dimensional pattern (see Fig. Ilip . 

We start with a generic situation corresponding to equation p3p or (j34l) . Since we are going to consider only small 
modulations of the soliton strings, we can follow the procedure suggested in 44] to get the Kadomtsev - Petviashvili 
equation for two-dimensional systems. 

Firstly we recall that equations (|3ip - (l32l) arc written in the moving coordinate systems. For the left- moving part 
of the solution in the laboratory frame of reference we have 



rr = VMJi- - Jz) \/32(l - Mo) d - 6^0 VMJi- - Jz) r 



1 T-i 



Mo 



J± - Jz \""^ V6 1 - Mo 
In writing the last equation we preserved the restriction = const 



g Jz ) rxxx 



(44) 
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The first term in the right-hand part plays the main role in the evolution of r^{X,T) and other terms represent 



small corrections with respect to the main contribution. According to [4J| we only need to calculate corrections to 
the main term coming from the slow modulation of the solution in the y-direction. This procedure gives us stable or 
unstable variants of the Kadomtsev - Petviashvili equation. The main term in the right-hand part of (j44p originates 
from the linear system ([T^]) of ([20)1 which can be easily written in the two-dimensional form by adding additional 



derivatives in the y-direction. What we need here is correction to the dispersion law ^ /c^ which can be written 



as 



kl + kl ^ ^J^{J^- J,) 32(1- f4) (^kx + 



for the left-moving part in our situation |105,]. As a result, the small modulations in the F-direction of solutions of 
(|44| can be described by the equation 



I T 



+5\/.'±(j±-A)v'32(i-f§)'-; 



YY 



rrx = \ VJ^iJi- - Jz) V^32(l - ^ig) ^i,^ „ 6^0 ^J^iJ^-Jz) ^x) x + 



in the moving coordinate system. 

Equation P5)) is the Kadomtsev - Petviashvili (KP) equation which describes the small transverse modulations of 
solutions of the KdV equation considered in the two-dimensional case. The stable Kadomtsev - Petviashvili equation 
corresponds to the same signs of the coefficients for Txxxx ^^'^ ''"yy ■ ^^is case the small modulation of a soliton 
string causes just the weak oscillations along the string and does not produce any instability. The opposite situation 
with different signs of the coefficients before r^xxxx ^-"^^ ^yy corresponds to the unstable situation where the soliton 
strings are unstable with respect to modulation along the y-axis. 

We can see then that the stable soliton string in two dimensions arises for the situation of equation ([55)) . i.e. 



./± - IJz J± - 7Jz 



which corresponds to the small values of Jz and the density n ~ 1/2 in the pattern. 
The solutions we considered in the opposite situation 



Jz>Jj7 or iMol > J^^^^^ 



are unstable from the point of view of the two-dimensional modulations. [106 

Let us say now that the analogous considerations can be performed also in the case of equation ([5^ so the results 
formulated above can be used also in the limit /iq — > 0. 
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We must certainly say that the Kad omts ev - Petviashvih equation is an integrable system from the point of view 
of the inverse scattering methods ( 2^ 10(^ ). The theory of equation (|45)) is very deep and brought many beautiful 



ideas in the theory of solitons. Let us just mention here two nice classes of solutions of in the stable and the 
unstable situation. 

1) The most interesting solutions of the Kadomtsev - Petviashvih equation in the stable situation are the two- 
dimensional A^-soliton solutions which are described in general by the formula 



U{X,Y,T) = ^{uj^T + k]^X + kl.Y + c^,...,uj^T + k'^X + k^Y + c'^) 

with some special functions ^{9^, . . . , 9^) ([so!]). 

The A^-solution solutions of the KP equation represent N plane interacting waves propagating at some angles with 
respect to each other. The interaction of the waves results in the phase shifts which can be rather big in the resonant 
case (@|)- 

2) For the unstable variant of the KP equation very interesting rational localized solutions ("lumps") can arise. The 
"lumps" represent localized both in X~ and F— direction solitons with rational dependence of coordinates. The 
interaction of solitons does not produce any phase shifts in this situation, so the solitons completely "forget" about 
each other after the interaction 

Let us emphasize here that the relation J± > was assumed everywhere in our considerations above and the 
properties we consider will be completely changed for the opposite situation J± < J^. Thus, as we pointed out 
already, the hydrodynamic approximation (jisp reveals an elliptic instability for the small values of fc (fc < 7r/2) in this 
situation which corresponds to a modulation instability of long-wave solutions of ([8]) in this case. In the same way, 
equation ()42[) becomes the focusing nonlinear Shrodinger equation in this situation which corresponds to the unstable 
behavior of the long-wave solutions of ([8]) either. However, the integrable nature of the focusing nonlinear Shrodinger 
equation leads to very interesting behavior of solutions also in this case. The most interesting part is the presence 
of the A^-soliton solutions for the focusing NLS equation which should be observed for > Ji. The corresponding 
two-dimensional equation for (|42p can be written in the form 

ihi^T = ^ {J^- J,) ijj - A{J^- J,)\iI^\'^iIj + Ah^ J^ij^xx + J^i'YY (46) 

The one-dimensional solutions of (|46p. however, are unstable with respect to the weak transverse modulations (Ji]) 
for Jz > J±. 



CONCLUDING REMARKS 



Soliton solutions in quantum systems is a subject of considerable theoretical interest. However, most of the earlier 
work focused on one dimensional systems, where special analytical tools, such as the Bethe ansatz solution, are 
available. For example, exact solitonic solutions were considered recently in a different quantum system in a series 
of papers [8-10]. Their analysis relied on the quantum inverse scattering methods, which are special to Id integrable 
systems. Our analysis in this paper is on constructing semiclassical solitons in two and three dimensional systems. 

States described by the wavefunction ([5]) correspond to collective excitations in the superfluid state. In the super- 
fluid state the U(l) symmetry is spontaneously broken, so the number of particles is not a good quantum number. 
Solitons which we discuss in this papers are semiclassical collective excitations. They can be thought of as spatially 
inhomogeneous coherent states representing non-linear excitations of the Hamiltonian. These solitons do not have a 
well defined number of particles. Within our approximations solitons have infinite lifetime. We expect that including 
coupling to other excitations may give rise to small but finite decay rate for the solitons, which may lead to dissipa- 
tive terms in the semiclassical dynamics. We expect that this should not change our conclusions qualitatively, since 
solitons should be robust against small dissipation (69| . 
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APPENDICES 

General approach for analyzing solitonic solutions in KdV-type equations 

The famous procedure of integration of the KdV equation ([is"]) is based on the connection of the KdV with 
the linear Shrodinger operator passing through the iso-spectral deformations according to the KdV evolution. The 
corresponding linear problems have the form 

- to = (47) 

for equation ([551) . and 

- i^xx - Utp = Ei: (48) 

for equation (|34l) . The connection of the KdV equations with the linear problems (|T7)) - (|^5)) gives a possibility to 
represent also equations ([55]) - (|34p in the equivalent form (|58|): 

±L ^ LA ~ AL (49) 



dT 



where the operators L, A have the form 



for equation and 



for equation (|34p . Representation p9|) of the KdV equation permits to consider the KdV evolution as the isospcctral 
deformation of the operator L using the exponent of the operator A as the corresponding basis transformation. 

According to the procedure represented in [2^ the scattering problem for the linear equations (|T7)) and P5|) plays 
the basic role in solving equations and in the rapidly decreasing case |C/(Ar)| 0, X ^ ±oo. Thus, if we 
consider the eigen-functions of ([T7)) or (HS]) having the asymptotic form 

^P{X) ~ e'''^ + b{k) e-'^^ , X ^ ^oo , ^P{X) ~ a{k) e^^^ , X ^ oo 

{k^ — E) and introduce the reflection and transition coefficients r{k), t{k) in the standard way we will have very 
simple evolution of the functions r{k,T), t{k,T): 

t{k,T) = t{k,0) , r{k,T) = e****'''^ r(fc, 0) 
according to the KdV evolution of U{X,T).^Q^ 
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In the same way, if the potential U{X) has bounded states '0„(X) with the energies En we will have En = const 
during all the KdV evolution. From the other hand, provided that the functions ^jJn{X) are normalized in the following 
way 

^„(X) ~ e'="^ , X^~oo , MX) ^ C„e-'="^ , X ^ oo 

{-kl = En) the evolution of the values C„(r) is given by C„(T) = e^^'^-'^ C„(0). 
The full set of the scattering data 

{r{k),En,Cn} 

gives the full information about the potential U{X) (0, 0, Is^l) such that the solution U{X,T) can be reconstructed 
at every time T using the values of r(fc, T), En, Cn(T). 

The potentials U (X) having zero reflection coefBcient r{k) = are called the reflectionless potentials and correspond 
to the exact A^-soliton solutions of the KdV-equation. The number of the bounded states (n = 1, . . . , A^) is equal to 
the number of solitons in the A^-soliton solution, so we can say that every bounded state in potential U (X) corresponds 
to a soliton in the solution U{X,T). The one-soliton solutions of the KdV-equation have the form 

2a2 

U{X,T) - ~ ^jj2^^^^4^3y^^^^ (50) 

for equation p3p and 



2a^ 

U{X,T) = (51) 

^ ' ch2(aX-4a3r + co) 

for equation (|34l) . Potentials (|50l) and (|5ip have exactly one bounded state according to linear problems (|T7)) and pS)) 
respectively with energy Ei — Ei{a) depending on the amplitude of a soliton. 

Analysis of solitons close to half-filling. Modified KdV equation 

In this section we discuss soliton solutions of two types of the mKdV equations : 

Ut + QU^Ux - Uxxx - (52) 
Ut + QU^Ux + Uxxx - (53) 

(we put a = here). 

Equation (1531) has two varieties of one-soliton solutions of arbitrary amplitude defined by the analytic formula 

± [^^==X + C 

Here U should be taken from one of the regions in the U-space where the value of expression vU'^ — is positive 
(see Fig. [IID. 

It is then easy to see that we can have either the particle-type or hole-type solitons, both moving to the right 
(v > 0) [lOSj l and connected by the transformation U ~U (Fig. [T3|) . 

Soliton velocity is proportional to the square of the amplitude v ^ and we can have arbitrary positive value of 
A. Explicit formula for the one-soliton solutions of (|53|) can be written in the form 

n 

u ^ ± 



ch {aX -a^T + ca) 
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FIG. 12: The left and right paths of integration w.r.t. U corresponding to hole- and particle-type one-sohton solutions of (l53l 




FIG. 13: The particle- and the hole-type solitons for equation (|53|) . 

Equation ([55)1 also admits more general soliton solutions. One can construct soliton solutions on a "pedestal" 
These solutions are defined by a more general analytic formula 



± 



dU 



X + c 



y/vU^ - [/4 - 2vUoU + 4[/3f/ + vU§ - 3U^ 

where two different paths of integration w.r.t. U are shown at Fig. [M) 

Again we can have solitons of the particle and hole type, both on a "pedestal" U = Ua moving with the speed v 
(our discussion is done in the moving frame) which can be represented by the following explicit formulas 



U = Uo + 



+ a2 ch (aX - {Wla + a^)T + cq) + 2[/o 



yjml + a2 ch (aX ~ [W^a + a^)T + cq) - 2Uo 



We have here v = 6C/q + while the amplitudes of the particle-type and the hole-type solitons are given by the 
formulas 



Ap.t. = 



V4C7JT^+2C/o 



^4{/2 + ^2 - 2C/o 
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FIG. 14: The left and right paths of integration w.r.t. U corresponding to the hole- and t particle-type soliton solutions on a 
"pedestal' Uo for 
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15: Solitons of the particle- and hole- type on the "pedestal" U = Uo > moving with the same velocity v for equation 



(see Fig. fTS]) . 

We can now see the difference in tlie particle- and liole-type solitons in this new situation. For C/q > the amplitude 
of a particle type soliton can be arbitrarily small for a — > 0, while the amplitude of the hole- type soliton is bounded 
from below by the value AUq (the situation is opposite for C/g < 0). We also see that solutions, which we consider, 
can be described as ordinary solitons of equation 



Ut + (6 C/q' + 12f/of/ + 6t/2) Ux + Uxxx = 

after the shift U U — Uq. This coincides with the general mKdV equation (jJH) after a Galilean transformation. 

We can claim then that regimes described by equation (p4|) (i.e. J± < 7Jz, or /Xq > {J± — 7Jz)/7{J± — Jz) if 
Jl > 7Jz) admit hole-type solitons after including the next nonlinear corrections. However, the small amplitude limit 
A is possible only for /io — for the hole-type solutions. As a result, we expect that new solutions, which we 
discussed above, can only be observed when 



J_L < 7Jz , Mo ^ 

and where changing from p4p to (j53p is quite natural. In the regime 




FIG. 16: The path of integration w.r.t. U corresponding to a one-soliton solution on a "pedestal" for equation (I52|) . 



fil > {J± - 7J,)/7{Ji_ - J.) 

it is easy to see that the limit /io ^ is possible only for Jl ^ 7J^. However, as we pointed out already, this situation 
is more complicated and should not be considered from the point of view of equations p4p or (j53p . Thus, we can see 
that hole-type solitons can arise in the regimes corresponding to equation p4[) for the situation Jl < 7Jz in the limit 
fiQ — T' +0 as a "reminiscent" of the region /io < as follows from the higher corrections to p4|) . 

The A^-soliton solutions of equation (|53|) as well as the solution of the initial value problem can be constructed in 
the form analogous to the case of KdV (see 64, 95]). 



We can see then that equation (|53|) gives a good limiting case of equation (p4)) for /ip — ^ in the situation Jl < 7Jz ■ 
Moreover, equation ((53)) provides a good limit for both cases hq > and fiQ < 0. The most remarkable feature of this 
regime is that both particle- and hole- type solitons with small amplitudes can coexist. The cubic nonlinear correction 
preserves the property of integrability of the corresponding evolution. Hence we expect that our analysis is applicable 
in the vicinity of the point /io = 0. 



Let us turn now to the regimes described by equation ([551) (i.e. < {J± — 7Jz)/7{J± — Jz), J± > 7Jz) which 
correspond to equation ((52|l for — Q. 

It is not difficult to see that equation ([5^ does not have real soliton solutions in ordinary sense and only the soliton 
solutions on "pedestal" can exist in this case. The one-soliton solutions on "pedestal" are defined by the analytic 
formula 



J ^-vU^ + U'^ + 2vUqU - AUi^U - vU^ + 

where the path of integration w.r.t. U is shown at Fig. 1161 

Explicit formula for the soliton solution can be written in the form 



f/ = ± 



2a^ 



v/C/q - ch {2aX - (UU^a - 8a^)T + cq) + Uo 



such that the soliton is of the hole- type for the positive " pedestal" and is of the particle- type for the negative pedestal 

(Uo > 0) (Fig. m- 

The amplitude of soliton 



26 





u 








Uo 








/ \ V 




/ 




X 




-Uo 



FIG. 17: The solitons of the hole type and of the particle type on the positive and negative "pedestals" for equation (|52p . 
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FIG. 18: The limiting form of a soliton solution for equation (|52|l . 



does not exceed the value 2Uo and can be arbitrarily small for a — > 0. The inverse scattering method and construction 
of the A^-soliton solutions on "pedestal" for equation ([52l) were considered in [79] and equation (|52p demonstrates 
that integrable properties are analogous to those of the KdV equation. 

We can see then that equation (|5^ gives a satisfactory limit of the regimes described by equation ([55)) (/Iq < 
(Jl — 7Jz)/7{J± — Jz), J± > 7Jz) in the limit iiq 0. We have to note, however, that the amplitude of solitons is 
restricted now by the value 2/io for /xq — > and soliton solutions disappear for fiQ = 0. Thus, generation of solitons in 
the regimes corresponding to equations ([55|) . ([5^ should be suppressed in the limit /Xq — J' 0. This should be contrasted 
to the regimes corresponding to equations ([M]). ([55)) . 

Appendix. Step decay close to half-filling 

One can use the inverse scattering method to solve initial value problems with localized initial perturbations for 
equations ([5^ or (155)) very similarly to what we discussed for equations (1551) or ([5^ . However, localized initial 
perturbation {U{X) 0, X ±00) will be a source of solitons at final stages only for equation ()53p for /ip = 0. 
The soliton part will be absent in the solutions of ([52)) . We also point out that for small /ip 7^ and big amplitude of 
initial perturbation for equation ([52)) (Vq >> ^iq) the "limiting" soliton (Fig. ITS)) in the limit T ^ 00 can arise (j75|). 

We also discuss briefly dynamics starting from the step-like initial state for equations ([52)) and ([53)) and the 
asymptotes of the corresponding solutions for T — > 00. According to the type of the solutions we considered above 
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FIG. 19: Appearance of oscillations for two different kinds of steps for equation (|39p in the case of different signs of U\ and 'U2- 
we will consider now the initial data such that 

\J{X) Ui , X ^ -00 , U{X) U2 , X ^ +00 

where both Ui and U2 are supposed to be small. 

Let us note first of all that the situation here is not pretty much different from those shown at Fig. [8] and Fig. |9]in 
the case when Ui and U2 have the same signs (say Ui, U2 > 0). So, the new features will arise here only in the case 
of different signs of Ui and U2 both for equations ([52]) and (|53l). 

Let us start again with equation (|53p . 

We have to say first that the oscillation region arises now for the both kinds of steps for the different signs of Ui and 
U2 (see Fig. [TO]) and the situation with just a decreasing of the steepness of initial data shown at Fig. [7] is impossible 
in this case. 

Both the situations shown at Fig. [19] for ([53| result in the generation of solitons on the final stage which have the 
particle type in the first and the hole type in the second situation (Fig. \20\i . 

We can see that the regimes of decay of step-like initial data for (|53|) include both the regimes coming from /xq > 
and fiQ < which is rather natural and gives a good limit for /ip ^ 0. 

Let us consider now the situation of equation (j52p corresponding to the small values of Jz and fiQ. Let us consider 
the initial data shown at the top of Fig. [T9l and suppose first that \Ui\ > \U2\. At the situation we describe the final 
stage of the oscillations development looks rather similar to that shown at Fig. |8] which is rather natural for the limit 
/io — ^ in the pattern. However, the limit \U2\ — ^ \Ui \ demonstrates quite new features here which are connected with 
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FIG. 20: The particle-type and hole-type solitons on the pedestals, Ui and U2, arising for two types of the step-like initial data 
for equation (I53p . 



the arising of a new solution for equation (f52|) . Indeed, for IC/2I \Ui\ the sohtons arising in the decay of the step-like 
initial data have a "limiting" form (Fig. [T8|) which is connected with the separation of the "shock- wave" solution 



U = -ath{aX ~ 2a^T + cq) (54) 

for |f/2| = \Ui\=a. 

Solution (|54|) plays an important role in the decay of the step-like initial data we consider for (|52|) for |J72| > 
Let us say that for general initial data having the form 



C/(-oo) = -C/(+oo) = a 

all the parts including ([5^ . solitons and the "wave-train" will generically arise ([Z^). 

It 's not difficult to understand also that for | C/2 1 > \Ui \ an additional step of the height | ?72 1 — | C/i | with the decreasing 
steepness will arise near the level U2 after the separation of solution (Fig. [^T|) . 

We have to say now that the step-like initial conditions of the second type (the bottom of Fig. [TO)) can be investigated 
just by the change U ^ —U. 

Let us mention here also the very interesting solutions of ([52]) including the soliton part and solution ([54)) . The 
soliton solutions coexist with solution ([M)) and the interaction of a soliton with ((54|) results in the phase shift and the 
soliton "flip" (Fig. [22]). 

Finally, we point out again that while considerations above were given for the function [/(X), representing Riemann 
invariants r^^'^^ {X), we can express the results in terms of physical variables p and k using equation (j24p (we also 
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FIG. 21: Additional step with decreasing steepness (T — >■ 00) arising after the separation of solution (|54[) for \U2\ > \Ui\. 

remind the readers that our analysis assumes the limit p ^ and fc — > 0). We find that for |?7i| > \U2\ (Fig. [531) the 
case 



corresponds to the case p{X) > 0, k{X) > 0. For \Ui\ — \U2\ ^ \Ui\ we also have p2 <^ pi, ki <^ k2 (see Fig. [23]). It 
is not difficult to see that conditions of this type can arise naturally after separating the right- and left-moving parts 
of initial conditions, as shown in Fig. [511 . 
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